function [ y ] = pis( theta, id)
load('mainfile1.mat');
load('matlab2.mat');
load('combnz.mat');

alpha=-0.1*theta(1);
c=theta(2);
%beta=theta(3:6); % <1x4 double>
mu=theta(7:31)*10;  % <1x25 double>
rho=theta(32:35);
%lambda=theta(36);
%gamma=theta(37);


n=find(consumer_id==id);
k=dimension(n(1));
m1=n(1);
m2=n(1)+k-1;
s= brandid(m1:m2)';

Si=matlab2(id,:); % matlab2=[6 21 9 0 0 0 0]; %i=1
tot=sum(logical(Si>0));

J=Si(find(Si>0)); % Si=[6 21 9];
Matrix=zeros(tot,1);

for h=1:1:tot
    Matrix(h)=(mu(J(h))+  x(n(1),:)*beta'+alpha*store_totprice(J(h)));    
end 
    %const1=-max(0,max(Matrix));
    %Matrix=Matrix+const1;
    Matrix=exp(Matrix);   
    powerset=cell2mat(combnz(1,tot-1));
    for l=1:(2^tot-1)
        temp(l)=powerset(l,:)*(Matrix);
        cost(l)=(sum(powerset(l,:)))*(c+x(n(1),:)*rho');
        powerset(l,:);
    end 
    
    temp(2^tot)=exp(0); %%%%%0+const1
    cost(2^tot)=0;
    total=log(1+(temp))-cost;
    %const2=-max(total);
    %total=total+const2; %% to avoide exp too large
    
sum2=0;
for i=1:1:k
    j=s(i);
    delta= mu(j) + x(n(1),:)*beta'+alpha*store_totprice(j);
    sum2=sum2+exp(delta);%%%summ=summ+exp(delta+const1);
end 
cost=k*[c+x(n(1),:)* rho']; %cost must be positive
mis2=log(1+sum2)-cost;
% mis2=mis2+const2;

y=exp(mis2)/sum(exp(total));
end

